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Summary 

A cell-vertex scheme for the Navier-Stokes equations, which is based on central difference 
approximations and Runge-Kutta time stepping, is described. Using local time stepping, 
implicit residual smoothing, a multigrid method, and carefully controlled artificial dissipative 
terms, very good convergence rates are obtained for a wide range of two- and three- 
dimensional flows over airfoils and wings. The accuracy of the code is examined by grid 
refinement studies and comparison with experimental data. For an accurate prediction of tur- 
bulent flows with strong separations, a modified version of the nonequilibrium turbulence 
model of Johnson and King is introduced, which is well suited for an implementation into 
three-dimensional Navier-Stokes codes. It is shown that the solutions for three-dimensional 
flows with strong separations can be dramatically improved, when a nonequilibrium model of 
turbulence is used. 

Introduction 

In recent years, considerable advancement has been achieved in the numerical solution of 
the Euler and Navier-Stokes Equations. Nevertheless, existing computer codes for the solution 
of the three-dimensional Navier-Stokes equations require large computing times. In order to 
make Navier-Stokes solutions useful in the flight vehicle design process, substantial improve- 
ments in efficiency and accuracy of the algorithms have to be made. 

Finite-volume methods based on explicit Runge-Kutta time stepping schemes have been 
shown to be very efficient in the solution of the Euler equations governing inviscid flows 
[1,2]. They are in widespread use now. More recently, they have been extended for a solution 
of the mass-averaged Navier-Stokes equations. The convergence of the basic scheme to the 
desired steady state solution, however, slows down considerably because of the time step lim- 
itation associated with the small mesh cells which are necessary to resolve the thin shear 
layers. This drawback of the explicit scheme can be overcome by applying several accelera- 
tion techniques, namely local time stepping, implicit residual averaging and a multigrid algo- 
rithm. Successful applications of these techniques have been reported in [3,4] for two- 
dimensional flows and there has been a partial success in three dimensions also [5]. 

In the present work a recently developed efficient and robust finite-volume multigrid 
scheme for the two- and three-dimensional Navier-Stokes equation is presented. The scheme 
employs a cell-vertex discretization rather than the usual cell-centered discretization because 
an analysis of the discretization errors indicates lower errors for meshes with high stretching 
or discontinuities in slope. The accuracy of the treatment of the far-field boundary is 
improved by representing the lifting surface by a vortex and using the induced flow field to 
determine the free-stream conditions. As the scheme is based on central differencing, an 
artificial viscosity model is used, which is particularly designed for high aspect-ratio cells. The 
boundary treatment of the dissipative terms, which is important for an accurate prediction of 
the skin friction, is discussed. A multistage Runge-Kutta time-stepping algorithm is used to 
advance the solution in time. The acceleration techniques, which are applied to obtain faster 
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convergence, are discussed. For an accurate prediction of turbulent flows with strong separa- 
tions, a modified version of the nonequilibrium turbulence model of Johnson and King [6] is 
introduced, which is well suited for an implementation into three-dimensional Navier-Stokes 
codes. 

Numerical results for subsonic, transonic and supersonic flows around airfoils and finite 
span wings are given. The accuracy and the convergence behavior of the scheme are investi- 
gated by varying the grid density and the coefficient of the artificial dissipation. The influence 
of the turbulence model on the numerical results is shown for several flow cases. It is 
demonstrated that the solutions converge very rapidly to the steady state. Typically, engineer- 
ing accuracy is obtained within 40 multigrid cycles for high Reynolds number transonic airfoil 
flows on a 321x65 mesh. For finite- span wings, converged solutions are obtained within 50-80 
iterations on a mesh with 288x65x49 points. 
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half span 
speed of sound 
skin-friction coefficient 
pressure coefficient 
drag coefficient 
lift coefficient 

Courant number of time stepping scheme 
Distance from the wall 
total internal energy 

unit vectors of Cartesian coordinate system 

coefficient of the second difference dissipation 

coefficient of the fourth difference dissipation 

Mach number 

unit vector of outer normal 

pressure 

Prandtl number 

cell face vectors in the T)-, and £- directions 
time 

temperature 

Cartesian velocity components 
volume 

vector of conserved quantities = (p pu pv pw pE) T 
Cartesian coordinates 
angle of attack 

ratio of specific heats, intermittency factor 
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Xfy X^, X ^ 

P 

a 

x 

CO 

Q 


central difference operator, boundary layer thickness 
coefficients of implicit smoothing according to (26) 
v. Karman’s constant = 0.4 

spectral radii of the Jacobian matrices associated with the t,-, T|-, ^-directions 

viscosity according to (4) 

density 

nonequilibrium factor in J.-K.-mod. turbulence model 
shear stress 

exponent in (12) and (27) 
magnitude of vorticity 
curvilinear coordinates 


Subscripts 

ij,k 

m 

max 

t 


w 


oo 


discrete quantity 

location of maximum of function g, equation (36), across boundary layer 
location of maximum of function F, equation (34), across boundary layer 
turbulent 
value at the wall 

value infinitely far away from the aerodynamic surface 


Governing Equations 

The integral form of the mass-averaged Navier-Stokes equations using nondimensional 
variables can be written as 

|jJJ*«V + ;jFlfdS=0 (1) 

V 3V 

where 

= (p pu pv pw pE) T . 

is the vector of conserved quantities with p, u, v, w, E denoting the density, the Cartesian 
velocity components and the specific total internal energy. The quantity V denotes an arbi- 
trary control volume with the boundary 3V and the outer normal if. The flux tensor F may be 
divided into its convective part F c and its viscous part F v as 


F = F C -F V 


F c = 


pu 1 X + pv i y + pw i z 
(pu 2 +p)!^ + puv ~i^ + puw 

puvT x + (pv 2 +p)l^, + pvwl^ 
puwl^ + pvw!^, + (pw 2 +p)T^ 
(puE+upJI^ + (pvE+vpji^ + (pwE+wp)T z 
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The equation of state for an ideal gas is used to calculate the pressure, p, and the temperature, 
T. That is, 

p = CH)p( E - (u 2 +v 2 )/2 ) , T = p/p , (2) 

where y denotes the ratio of specific heats. The elements of the shear-stress tensor and the 
heat-flux vector are given by the constitutive equations for a Newtonian fluid as follows: 

a xx = 2 P u x " 2/3|i(u x +v y +w z) 

CTyy 2p.Vy 2/3 P(U x +Vy+X Z ) 

(T zz 2pv z 2/3p(u x +v y +x z) 

a xy = °yx = ^(V V x) (3) 

°xz = °zx = ^(U z +W x ) 


Oy Z — 0 Z y — P(v z +Wy) 

q x = -k5T/3x, q y = -k3T/9y, q z = -k9T/8z . 


The nondimen sional viscosity, |l, is assumed to follow an empirical power law 


Y 1/2 M„ 

Re„ 


P = 


(T/Tj 


, 0.75 


( 4 ) 


and the heat conductivity is 

k = _)L - \L . 

Y~1 Pr 


( 5 ) 


For turbulent flows, the laminar viscosity, p, is replaced by p+p t and p/Pr is replaced by 
p/Pr+p/Pq, where the eddy viscosity, p t , and the turbulent Prandtl number, Pr t , are provided 
by a turbulence model. The turbulence models used in the present work are described in a 
subsequent paragraph. 


Spatial Discretization 

The discretization of (1) follows the method of lines, i.e. discretization in space and time 
are done separately. The domain around the aerodynamic body is partitioned with hex- 
ahedrons. The discrete values of the flow quantities are located at the vertices of the mesh 
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cells. The eight cells surrounding a vertex form a super cell as shown in Fig. 1. The integral 
equation (1) is approximated by the spatial discretization yielding 



In the following, the net inflow of mass, momentum and energy associated with the convec- 
tive terms (^ c j j k , the viscous terms (^ v y k and the artificial dissipative terms Dy * are dis- 
cussed. For this purpose the super cell around the point (i,j,k) for the simplified case of a 
plane mesh is shown in Fig. 2. 

The surface integral of (1) is evaluated for each component cell using an arithmetic aver- 
age of the flux quantities at the vertices to determine the values on each of the boundary seg- 
ments. Then the resultant convective inflow of mass, momentum and energy associated with 
point (i j,k) is computed by summing the contributions of the component cells. According to 
[7], the scheme is at least first-order accurate, if two conditions are satisfied. (1) The distribu- 
tion of the normal vector over each surface segment has to be smooth, that is, if the grid is 
refined, the four grid points defining a segment tend to lie within a plane. (2) The shape of 
each surface segment has to approach a parallelogram with grid refinement. Note, that these 
conditions still allow discontinuities in the slope of the grid lines and discontinuities in the 
rate of grid stretching. On completely smooth meshes, the discretization is second-order accu- 
rate. 


Next, the viscous fluxes required to determine the solution at the point (i,j,k) are approxi- 
mated using the auxiliary cell with the dashed boundary shown in Fig. 2. The viscous fluxes 
contain first derivatives of the flow variables, which are computed using a local transformation 
from Cartesian coordinates to the curvilinear coordinates 5 ,tj,£, i. e. 

= 9<t> 35 9<j>dr| d£ d£ 

Vx as a* dn dx ac dx 


( 7 ) 


The derivatives <j>^, ^ and <j>^ are approximated using finite differences, whereas the cell face 
vectors 5^, S^, 3?^, which are written as 

^ = (S^ x S^, S^) T 

= ($Tix S^ y 

^ = c S^) T , 


and the volume V are used to compute metric derivatives. If <)> x is to be approximated at 
(i+l/2j,k) we obtain 


VHVi+l/2,j,k = 77 (8) 

v i+l/2.j.k 


where 5^, 5^, 6^ denote central difference operators in the curvilinear coordinate directions. 
In practice, coordinate grids for the resolution of viscous shear layers are highly stretched in 
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the direction normal to the layer and therefore, an one-dimensional error analysis of the 
discrete approximation to the viscous terms can give useful insight. In [4], it was shown that 
the present scheme is first-order accurate on general stretched meshes and second-order accu- 
racy is obtained on smoothly stretched meshes. In the present version of the code, the viscous 
terms have been simplified by taking into account gradients in the direction normal to the 
viscous shear layers only (thin layer approximation). 

In order to prevent odd-even point decoupling and oscillations near shock waves, 
artificial dissipative terms are added to the governing discrete equations. The artificial dissi- 
pation model considered in this paper is based on the work of Jameson, Schmidt and Turkel 
[1]. A blend of fourth and second differences is used to provide third-order background dissi- 
pation in smooth regions of the flow and first-order dissipation at shock waves, and is given 
by 

By* = (D 2 5 - D 4 5 + D 2 , - D 4 , + D 2 ? - D\) vjy ik . (9) 

The second and fourth difference operators read 

D 2 4 tfiO* = 'VV i + l,2J.k'e <2, i.l,20*)^V>iJJc . (10) 

D 4 . w/2J.k' e<4) i*i'20.k) A 5 V ^i.i* • (11) 

where A^, are forward and backward difference operators in the ^-direction. In order to 
avoid excessively large dissipation levels for cells with high aspect ratios and to maintain 
good damping properties of the scheme, a variable scaling factor of the dissipative terms is 
employed, which is a function of the spectral radii of the Jacobian matrices associated with 
the ^,r|,£ directions and accounts for varying cell aspect ratio 

^ i+l/2j,k = ^ i+l/2,j,k ' ^i+lttj.k > (12) 

^i+l^j.k — 1 "h m ^ x ((^r| i+l/2,j,k^i; i+l/2j,k) * i+l/2,j,k^^ i+l/2,j,k) ) * 

^ i+l/2,j,k — ^i+l/2j,k ' ^ i+l/2j,k ^ c ^ i+l/2,j,k ^ » 

K\ i+l/2.j,k = ll5> i+l/2,j,k * ^T| i+l/2,j,k 1 + C 1^ j + i/2,j,k 1 » 

\ i+l/2,j,k = l1[ ^i+l/2j,k ' ^ i+l/2,j,k 1 + C 1 ^ i + i/2.j.k 1 » 

where I?j + i/ 2 j,k denotes the velocity vector and c is the speed of sound. The use of the max- 
imum function in the definition of O is important for grids where and X^/X^ are very 

large and of same order of magnitude. In this case, if these ratios are summed rather than tak- ( 
ing the maximum, too large dissipative terms are obtained, which may destroy the conver- 
gence of the time-stepping scheme. In general, coordinate grids around three-dimensional 
configurations exhibit larger variations of the cell aspect ratio than two-dimensional grids. 
Therefore, we have observed somewhat less freedom in the choice of the exponent O) in three 
dimensions. It has been found that the choice <a=0.5 yields a robust three-dimensional scheme. 
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The coefficients e (2) and e (4) use the pressure as a sensor for shocks in the flow field. They are 
defined as 

^ \+i/2,j,k — ^ ) max(v i _ ]0ik ,Vi Ji k,v i+liJik ,v i+ 2 Ji k) (13) 

_ | Pi— 1 .j.k ~ 3pj,j,k + Pi+l,j,k | 

'' i,k Pi-lJ.k + 2 Pij,k + Pi+l.j,k 
e<4) i+i/2j.k - max< 0 , (k< 4 > - e (2 > w;2J ,k)> 

where k (2) and k^ 4) are constants. The dissipation operators in the r\- and ^-directions are 
defined in a similar manner. 

Boundary Conditions 

At subsonic inflow/outflow boundaries, locally one-dimensional flow normal to the boun- 
dary is assumed. According to [8], the governing flow equations are linearized around the 
values at the end of the previous time step, and the characteristic variables corresponding to 
outgoing characteristics are extrapolated from the interior. The characteristic variables 
corresponding to incoming characteristics are determined from the free stream. 

For two-dimensional flows around airfoils, it is well known that the free-stream condi- 
tions have to be specified by taking into account the circulation around the airfoil to obtain 
correct values of lift and drag for a given angle of attack. The influence of an improper 
determination of the free stream at the far field may be felt even if the distance to the far field 
is as large as 20-50 chord lengths. The usual practice in calculating the free-stream condi- 
tions for flows around finite-span wings has been to use simply the undisturbed onset flow 
conditions. In fact, if the far field is sufficiently far away from the wing surface, the distur- 
bances created by the wing behave like those coming from a point singularity rather than a 
line singularity as in two-dimensional flow. Hence, a disturbance would decay more rapidly 
in three-dimensional flow than in two-dimensional flow, thus allowing the use of a smaller 
distance to the far-field. 

More accurate results can be expected if the free-stream conditions used in three- 
dimensional flow computations take into account disturbances created by the aerodynamic sur- 
faces. Klunker [9] shows that the small-disturbance potential at the far field is dominated by 
terms representing the lift of the aerodynamic surface. He derives an equation in which the 
potential at the far field is determined by integrating the contributions of singularities 
representing the local circulation over the entire wing surface. This approach has been used in 
[10] to determine the free-stream conditions in an Euler code. 

If free-stream conditions are to be specified for general configurations a numerical 
integration over the aerodynamic surface for each point in the far field seems to be untract- 
able. In particular, if multiblock structured meshes are to be used, this integration requires 
considerable amounts of programming effort and storage. The determination of the free 
stream is much simpler if the wing is replaced by a horse-shoe vortex. The components of the 
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disturbance velocity vector induced by a horse-shoe vortex in compressible flow are* 


p 2 r 

u = - c — 


v = 


271 
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2TC X 2 +p 2 y 2 
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Vv+ Vvl 
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( 1 +- 2 -) 

(z+a) 2 +y 2 ' Vv+ 7 (z-a) 2 +y 2 Vvl 


■( 1 + 




(14) 

(15) 

(16) 


where 


[3 = Vi -Mj, 

y + = x 2 +P 2 (z+a) 2 + P 2 y 2 , 
\j/_ = x 2 +p 2 (z-a) 2 + P 2 y 2 . 


The circulation of the vortex is related to the lift of the wing by 



c C 

v mean v " 


L > 


(17) 


the quantity a denotes the half span of the horse-shoe vortex and c mean is the mean chord of 
the wing. The position of the horse-shoe vortex with respect to the Cartesian coordinates x,y,z 
is shown in Fig. 3. The induced velocity is infinite downstream of the wing where the vortex 
crosses the boundary. In reality, however, this is by no means a special point in the flow field. 
Furthermore, numerical experience indicates that the treatment of the outflow boundary down- 
stream of the wing has no significant influence on the results. Therefore, the following numer- 
ical treatment has been employed. The denominators in (15) and (16) which are responsible 
for the singular behavior of v and w at downstream boundaries, are limited in their value for 
x>0, i.e. (z+a) 2 +y 2 and (z-a) 2 +y 2 are kept at values larger than a/2 for x>0. This means that 
the induced velocities are artificially reduced in their value within the distance a/2 around the 
vortex line which crosses the downstream far-field boundary. 

For all grid points at the far field, the geometric relations in the formulas given above 
can be evaluated once and stored. At each time step, these influence coefficients are multi- 
plied by the total lift of the wing and added to the onset flow. Thus, the work to obtain the 
free-stream velocities is negligible. The pressure and the density of the free stream are then 
obtained using the assumption of isentropic flow and constant total enthalpy. The numerical 
results obtained with the present determination of the free-stream values show a reduced sen- 
sitivity to the far-field distance (see section on numerical results). 


* These relations have been derived in unpublished ’Notes on Linearized Subsonic Wing 
Theory’ by E.B. Klunkcr and K.C. Harder. 
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At solid walls, the no- slip condition is enforced by setting 

u = v = w = 0 (18) 

The continuity and energy equations are solved at the grid points lying on the surface assum- 
ing an adiabatic wall. For this purpose, a control volume is formed with the four nearest cells 
outside and their mirror images inside. When updating the convective terms, the flow vari- 
ables at the computational points inside are obtained as the symmetric images of the values 
outside. For the viscous terms, the density and the total energy inside are obtained as sym- 
metric images, whereas the velocity components are taken as the antisymmetric images of the 
values just outside. 

It has been found that the computed velocity distributions near the wall may be 
significantly affected by the artificial dissipative terms if the stencil of the dissipative terms 
normal to the wall is not properly defined. The reasons for this are the high grid stretching 
normal to the wall and the steep velocity gradients of the turbulent boundary layer. When 
solving the momentum equations at the grid point just outside the wall at j=3, the dissipation 
operator needs the dependent variables just inside. The extrapolation formula 

V>l,l,k = = 2 (19) 

rather than the usual linear form 

= A^iaic (20) 

results in much less sensitivity of the solution to artificial viscosity and has had no significant 
drawback on convergence in the computations made so far. 


Time Stepping Scheme 


The system of ordinary differential equations which is obtained by the discretization in 
space is advanced in time with a five-stage Runge-Kutta scheme. A hybrid scheme is used 
where the physical viscous terms are computed on the first stage and frozen for the remaining 
stages and the artificial dissipative terms are evaluated on the first, third and fifth stages of the 
scheme. At the (m+l)st stage, 


#<""•'> = $< 0 ) _ 


n=0 


( 21 ) 


^(m) = 3 c (m) + ^ (m) 


where V^ (0) is the solution at the old time level, and At is the time step. The subscripts c and 
v refer to convection and physical diffusion contributions to the discrete flow equations. The 
stage coefficients, a m are taken as 


a! = 1/4, 0 C 2 = 1/6, a 3 = 3/8, a 4 = 1/2, a 5 = 1 . 


( 22 ) 
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The weighting factors of the artificial dissipation satisfy the condition 

XYmn = 1 • (23) 

They are 

Yoo = 1 

Yio = 1 . Yn = 0 , 

Y20 = 0 -Y 3 ) , Y 21 = 0 , Y 22 = Y3 * (24) 

Y 30 = (!-Y 3 ) . Y31=0. Y32 = Y 3 » Y33 = 0» 

Y 40 = (I-Y 3 KI-Y 5 ) . Y41 = 0» Y42 = Y3(1-T5) . 

Y 43 = 0 . Y44 = Y5. 

where y 3 = 0.56, y 5 = 0.44. It has been shown in [3], that this scheme has a particularly large 
parabolic stability limit. The steady-state solution is independent of the time step and there- 
fore, the scheme is amenable to convergence acceleration techniques. 

Three methods are employed to accelerate convergence of the basic explicit time step- 
ping scheme. These techniques are as follows: (1) local time stepping, (2) implicit residual 
smoothing and (3) multigrid. 

With local time stepping, the solution at each mesh point is advanced at the maximum At 
allowed by stability. Both convection and diffusion stability limits are included in the compu- 
tation of At as 

V 2 

At = CFL — 5 - 5 -“ , (25) 

V(^ + ^n+^) + c di ffV diff <S^S^ + V-n + ^-S^) 

where 

. 4 y . 4+m 

v^maxCp -£)— , 

and the coefficient c diff = 2 is used. 

The stability range of the basic time stepping scheme can be extended using implicit 
smoothing of the residuals [11]. For three-dimensional flows the residual smoothing is applied 
in the form 

(l-e^Kl-e^V^XI-ejV^) R m = R„ (26) 

where R m is the explicit residual at the Runge-Kutta stage m in (21), and Rnj is the final resi- 
dual at stage m after the sequence of smoothing in the r\ and £ directions. 

The use of constant coefficients in the implicit treatment has proven to be satisfactory 
(extending the CFL number by a factor of two to three) even for meshes with cells of high 
aspect ratio, provided additional support such as enthalpy damping [1] is introduced. How- 
ever, the use of enthalpy damping, which assumes constant total enthalpy throughout the flow 
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field, precludes the solution of problems with heat transfer effects. The need for enthalpy 
damping can be eliminated by using variable coefficients e^, e^, and that account for the 
variation of the cell aspect ratio. 

For this purpose, consider a cell where the edge lengths in the and £- directions are 
much longer than that in the ri -direction. The explicit time step is limited by the characteris- 
tic wave speed in the direction of the short cell edge. It is obvious that, for an extension of 
stability, implicit smoothing is required in the T]-direction. If the same implicit residual 
smoothing is also applied in the other directions, where the characteristic wave speeds are 
much smaller than the stability limit, the damping behavior of the scheme, which is optimal at 
wave speeds near the stability limit, is impaired. To overcome this problem, Martinelli [3] 
has given a formula, in which the smoothing coefficient is a function of the characteristic 
wave speeds. Martinelli’ s form and variable smoothing coefficients in general are discussed in 
detail in the unpublished work of Swanson. Here, it is extended to three dimensions, i.e. 
may be defined as 


1 PFT Xr „ 

e*. = max(0 , 4( - t V — (l+maxCC-^-r^-^r))) 2 -!) 

^ 4 CFL explicit XjH-max^,^) X^ X^ 


(27) 


Finally, a multigrid method is employed which is based on the work of Jameson [2]. For 
the multigrid process, coarser meshes are obtained by eliminating every other mesh line in 
each coordinate direction. The solution is transferred to coarser meshes by injection. Residu- 
als are transferred from fine to coarse meshes by a weighted average over the fine mesh grid 
points which are nearest to the point on the coarse mesh. A forcing function is constructed so 
that the solution on a coarse mesh is driven by the residuals collected on the next finer mesh. 
This procedure is repeated on each succeeding coarser mesh until the coarsest mesh is 
reached. Then, the corrections are transferred to the next finer mesh by trilinear interpolation. 
A fixed W-type cycle is used to execute the multigrid strategy. The robustness of the mul- 
tigrid scheme is greatly improved by smoothing the total corrections before the solution is 
updated. That is, 

$(n+i) _ $ n +A$ tot (28) 

where 

A$ tot = A^ f + AV^ C . 

The quantity A$ F is the solution correction from the finest grid, and AV^ C is the resultant 
solution correction from the coarse grids. The smoothing reduces high frequency oscillations 
introduced by the linear interpolation of the coarse mesh corrections and hence, convergence 
of the scheme is obtained for a much broader range of dissipation coefficients. The factored 
scheme of (26) with constant coefficients (8^ = 8^ = = 0.2) is used for this smoothing. 

Additional robustness and efficiency of the multigrid scheme is achieved by computing more 
than a single time step on the coarse meshes. In the first place, the use of multiple time steps 
on the coarse meshes improves the convergence, because the solution is further advanced in 
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time on the coarse meshes. Secondly, multiple time stepping improves the damping charac- 
teristics on the coarse meshes. Thus, smaller values of coarse-mesh dissipation may be used 
which results in an improvement of the convergence to the steady state. Finally, the applica- 
tion of a Full Multigrid method provides a well conditioned starting solution for the finest 
mesh being considered. 

Turbulence Models 

In the present work two turbulence models are considered. The first model is the widely 
used algebraic model of Baldwin and Lomax [12], which is based on the eddy viscosity 
hypothesis. The eddy viscosity is defined as 

= min(|i t i , |i Uo ) . (29) 

In the inner region of the boundary layer, 

Mi,i = p(K D d) 2 Q (30) 

where k is von Karman’s constant, Q is the magnitude of the vorticity, and the coordinate d 
is the distance from the wall. The quantity D is the van Driest damping factor 

D = 1 - exp(-y + /26), y + = d . (31) 

M-w 

Note, that the maximum laminar shear stress across the layer, X lmax rather than the shear 
stress at the wall is used to define the damping term in order to prevent vanishing eddy 
viscosity at the separation points. In the outer region, 

|i t ,o - 1.6 • 0.0168pF wake y , (32) 

where 

F\vake = m i n (dmax^max ’dmaxUdif /Fmax) • (33) 

F max is the maximum value of 

F = D Q d (34) 

across the layer, and d max is the value of d at which it occurs. The quantity U dif is the 
difference of the flow speed across the shear layer. The intermittency factor y is 

Y= 7 (35) 

l+5.5(0.3-d/d max ) 6 

This turbulence model will be denoted by B.-L.. 

Frequently, it has been observed that in separated flows, the assumption of an equili- 
brium of the turbulence is not valid. Both convection and diffusion of turbulence have to be 
taken into account to allow for a more accurate determination of the turbulent stresses. For 
this purpose, Johnson and King [6] have devised an extended eddy viscosity model for 
separated wall boundary layers. They postulate that the viscosity of the inner layer near the 
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wall and the shape of the eddy-viscosity distribution in the outer layer may be accurately 
described with algebraic relations assuming equilibrium. The equilibrium distribution of the 
outer layer, however, is multiplied with a nonequilibrium factor in such a way that a transport 
equation for the maximum shear stress in the layer is fulfilled. Thus, the complete model is a 
nonequilibrium model of turbulence. In the following, a slightly modified version of this 
model is described. The model is particularly convenient for use in three-dimensional Navier- 
Stokes codes by using the Baldwin-Lomax formulation for the outer viscosity as proposed in 
[13]. 


First, instead of using the maximum turbulent shear stress as originally proposed in [6], a 
turbulence reference quantity. 


g m = (Q m/p) m 


(36) 


is introduced, where the index m denotes the maximum of g across the shear layer. Note, that 
g m is invariant to the coordinate system used, and therefore, it can be easily used in Navier- 
Stokes codes. According to [6], the eddy viscosity is defined by an exponential blending of 
the inner and the outer viscosity 


Pt = Pt,o 


1-exp 


Pt,i 

P*> . ‘ 


(37) 


The inner viscosity is written as 


Pt,i = pD 2 K d g m 1/2 , 


D = 1-exp 


-d Vpw max (p w g m , T w ) 


A p w 


(38) 


The use of the maximum function in (38) is important for numerical stability on coarse 
meshes, where the turbulent boundary layer is not well resolved. According to [14], a value of 
17 has been chosen for A + . In the outer layer, the original model of Johnson and King uses 
the flow velocity at the edge of the boundary layer and the displacement thickness to define 
the outer viscosity. The displacement thickness, however, is not easily determined in three- 
dimensional flows. In order to avoid numerical problems associated with searching through 
the boundary layer in three-dimensional flows, the viscosity in the outer layer is determined 
by the relation of Baldwin and Lomax 


Pt, o = a 1.6 • 0.0168pF wake Y . (39) 

The quantity a is the aforementioned nonequilibrium factor, and F wake is defined according to 
(33). The intermittency factor, however, is redefined as 

= 1 
Y l+5.5(d/5) 6 


( 40 ) 
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where 5 denotes the thickness of the boundary layer. In [15] it has been demonstrated that 


5=1-9 d max 


(41) 


is a good approximation for wall boundary layers, and this relation is incorporated into the 
model. Furthermore, the behavior of the model in separated flow can be improved, if the dis- 
tance d in (34) for the function F is computed according to Fig. 4 rather than taking the dis- 
tance from the wall. 


Assuming that g is proportional to the kinetic energy of turbulence, the distribution of g m 
over the aerodynamic surface is determined by a differential equation given in [6], For three- 
dimensional flows, the equation may be written as 


dgm 

dt 


d gm dg. 


m 


+u »ir +v ”-3r +w ” 3z 


3g m a l 8m 


__l — : — — (g g 

^ t VSm t>eq,m 


1 ' 2 )+gm 3 ' 2 D m = 0 


(42) 


Here, u m , v m , and w m are the Cartesian velocity components at the location of g m . The index 
eq,m denotes the equilibrium value at the location of g m , and the length scale, L^, is defined 
as 


L m = for c^/8 < 0.225 , (43) 

L m = 0.095 for c^/5 > 0.225 . 


The diffusion term, D, is defined as 


a 2 max(0, a 1/2 -l) 
5(0.7 - (d/S)J ‘ 


(44) 


The model constants apO.25 and a 2 =0.5 are used here. Equation (44) states that the diffusion 
of turbulence is neglected in retarded flows, where a<l, whereas in recovering regions, the 
diffusion is assumed to be proportional to a 1/2 -l . If g m is replaced with g m =g m _1/2 > a linear 
equation is obtained as 


dgm dg m 

3t +Um 9x +Vn 


dg : 


m 


dg, 


m 


l l 


ay dz 2 L_ 


1/2 a 


a i 


'geq.m gm ^ (D+ ^ ) 0 . 


(45) 




(45) is valid along a surface which is determined by the location of the maximum of the refer- 
ence quantity g. For numerical convenience, (45) is solved on the surface grid of the aero- 
dynamic body rather at the location of g m . The misalignment between these surfaces creates 
errors in the convective terms of (45). These errors are reduced by replacing the velocity com- 
ponents u m , v m , and w m in (45) by their projection onto the body surface. It is believed, that 
the remaining errors, which are due to the skewness of the grid lines crossing the shear layer, 
are small in comparison with the uncertainties present in the basic assumptions of this tur- 
bulence model. Note, that the time derivative in (45) is retained, so that a time stepping 
scheme can be used for the numerical solution. Using Gauss’s theorem, a cell centered 
finite-volume discretization with central differencing is applied to (45). Because g m does not 
vary normal to the surface, only a layer of one cell thickness around the surface is required to 
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solve (45). Fourth difference dissipative terms are used to damp the discrete equations. The 
explicit five-stage Runge-Kutta scheme (21) is used for time stepping. The linear source term 
in (45) is treated implicitly, and thus, CFL-numbers around 3 can be used. The convergence 
to the steady state is further enhanced by local time stepping. For the computation of flow 
cases with strong separation the turbulent viscosity should be rapidly adjusting to the flow 
field in order to allow for convergence. Therefore, 10 explicit time steps are executed at each 
time, when the turbulence model is updated. The computational time for the solution of (45) 
is very small in comparison with the total effort, because only a single equation for a quantity, 
which does not change normal to the surface, is involved. 

Once values of g m have been computed by the numerical solution of (45), they can be 
compared to values of g m from (36). The ratio of both these values is then used to update o, 
as described in [13]; and eventually, convergence is obtained. The nonequilibrium turbulence 
model described herein will be denoted with J.-K.-mod. . 

Numerical Results 

The cell-vertex multigrid scheme described here has been applied to a wide variety of 
two- and three-dimensional flows. The two-dimensional subsonic laminar flow around NACA 
0012 airfoil at the free-stream Mach number, = 0.5, the angle of attack, a = 0°, and the 
Reynolds number based on the chord length Re M = 5000 is considered first. A C-type coordi- 
nate mesh of 257x65 points is used, where 129 points are distributed over the airfoil surface. 
The first spacing away from the wall is 5X10 -4 chord lengths. Fig. 5 shows the convergence 
history, the streamlines, and the distributions of pressure and skin friction. The computed 
values of pressure drag and skin-friction drag and the separation point agree well with the 
results of Swanson and Turkel [16]. The solution converges rapidly, and, noting additional 
convergence criteria given in Table 1, a sufficiently accurate solution is obtained in less than 
40 multigrid cycles on the fine mesh. 

The transonic turbulent flow over the RAE 2822 airfoil is used to demonstrate the con- 
vergence behavior and the accuracy of the method for high Reynolds number turbulent flows. 
A mesh of 385x65 points is used, where 257 points are distributed over the airfoil surface, 
and the first spacing away from the wall is 10“ 5 chord lengths. Frequently, a mesh with 
321x65 points has been used also, which has the same distribution of grid points at the sur- 
face, but less points along the wake line. In Fig. 6, the convergence behavior is displayed for 
M„ = 0.73, a = 2.79°, Re^ = 6.5xl0 6 . Fig. 6a shows the influence of the fourth-difference 
dissipation on the convergence of the solution. The robustness of the scheme is demonstrated 
by the fact that the dissipation coefficient is varied by an order of magnitude without destroy- 
ing the convergence of the scheme. Note, that the results of Fig. 6a have been obtained using 
single time steps on the coarse meshes. The effect of multiple time steps on the coarse 
meshes is shown in Fig. 6b. If two time steps are used, the convergence is significantly 
improved, and a residual reduction of 10 orders of magnitude is obtained within 340 multigrid 
cycles on the fine mesh. The convergence criteria in Table 1 indicate that sufficiently 
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converged solutions are obtained within 40 multigrid cycles. Fig. 6c shows the convergence 
with the nonequilibrium model of turbulences, which is slightly slower than with the B.-L. 
model. 

Next, the accuracy of the scheme is examined using a variation of the grid density. For 
this purpose, a coarser mesh of 193x33 points is created by omitting every other point in each 
coordinate direction. Additionally, a very fine mesh of 577x97 points has been used. The 
variation of the coefficients for lift, pressure drag, and friction drag with number of mesh 
points N is presented in Fig. 7. The effect of the second and fourth difference dissipation is 
also indicated. On coarse meshes, the discretization error is obviously dominated by artificial 
dissipation. The finer meshes allow the extrapolation of the coefficients to their values for an 
infinitely fine mesh. For the 385x65 mesh (N -1 =4.0xl0 -5 ), the predicted lift is within 1.5 per- 
cent, the pressure drag is within 3 counts, and the friction drag is within 0.3 count of the 
extrapolated values. For the 577x97 mesh (N -1 =1.8xlO~ 5 ), the predicted lift is within 0.5 per- 
cent, the pressure drag is within 1 count, and the friction drag is within 0.1 count of the extra- 
polated values. Fig. 8 shows pressure and skin-friction distributions for different grid densi- 
ties, with a comparison to the experimental data of [17]. The main features of the flow are 
essentially captured with the 193x33 mesh. The differences between the solutions on the 
385x65 and 577x97 meshes are small. 

In Fig. 9, the effect of the turbulence model is presented. For both the equilibrium B.-L. 
model and the nonequilibrium J.-K.-mod. model, the position of the shock and its strength 
agree fairly well with the experimental data. The measured lift coefficient, however, is in 
better agreement with the results of the J.-K.-mod. model than with B.-L. The skin friction 
distribution of the J.-K.-mod. computation shows relatively high values downstream of the 
shock, which are not observed experimentally. The high skin friction in the recovery region is 
caused by the definition of the inner viscosity, (38), which is scaled by the maximum of the 
quantity g rather than taking local quantities as in (30). The definition (38) should be refined 
in future work. 

In the following, the behavior of the turbulence models is presented for test cases with 
strong separation, which have been used recently to compare various computer codes for 
viscous transonic airfoils [18]. In Fig. 10, the flow around the RAE 2822 airfoil at 

= 0.75, a = 2.81°, Re,*, = 6.5xl0 6 is considered. With the B.-L. model, the shock location 
is predicted 8% of the chord length downstream of the experimental data. The agreement of 
computation and experiment is considerably improved, when the J.-K.-mod. model is used, 
however, there are still some differences remaining. The bump in the pressure distribution 
downstream of the shock, as obtained with the J.-K.-mod. model, indicates that the height of 
the separated region, as well as the recovery process downstream are somewhat overpredicted. 

Similar trends are observed for the next test case, the NACA 0012 airfoil at = 0.799, 
a = 2.26°, Re^ = 9.0xl0 6 . Fig. 11 shows much better agreement in pressure distributions and 
lift and drag values, when the J.-K.-mod. model is used instead of the B.-L. 
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Drag polars have been computed for the NACA 0012 airfoil at M„ = 0.7 and 
Re*, = 9.0x1 0 6 . The results of these computations are shown in Fig. 12 (lift versus angle of 
attack) and Fig. 13 (lift versus drag). Experimental results of Harris [19] with wind-tunnel 
corrections included are also displayed. Good agreement of computations and measurements is 
obtained at lower angles of attack. At high angles, the B.-L. model yields considerably higher 
lift values, than the J.-K.-mod. In particular, the J.-K.-mod. model predicts maximum lift with 
Cl, max = 0-728 at a = 6°, whereas the maximum lift is C Lmax = 0.87 at a = 7° with the B.-L. 
model. There is no experimental value of the maximum lift for this particular Mach number 
given in [19]. The experimental values of the maximum lift at Mach numbers below and 
above, C L>max = 0.82 for M„ = 0.65 and C Lmax = 0.55 for M«, = 0.76 indicate that the result 
of the J.-K.-mod. model is close to the probable result of the wind tunnel. 

The last two-dimensional test case is the supersonic turbulent flow around NACA 0012 
airfoil at M m = 2.0, a = 0°, Re«, = 10.0xl0 6 . The convergence history and the distributions of 
pressure and skin friction are shown in Fig. 14. Generally, good convergence rates are 
obtained for Mach numbers up to 2. For higher Mach numbers, no solution could be obtained 
with the present scheme. In Fig. 15, the Mach contours and the mesh near the airfoil are 
shown. No attempt has been made to adapt the mesh to the bow shock, and therefore, the bow 
shock is not very well resolved in this computation. 

The three-dimensional version of the code has been applied to the flow over the finite- 
span ONERA-M6 wing and to the three-dimensional flow over an airfoil mounted in the wind 
tunnel, where the viscous layer along the wind-tunnel side wall has been included in the simu- 
lation. In the present report, only the computations of transonic flows over the ONERA-M6 
wing are presented. The results for the interaction between a wing and a side-wall boundary 
layer are given in a subsequent paper. 

The computational domain around the wing is discretized using C-type topology in the 
streamwise direction and an O-type topology in the spanwise direction with 289x65x49 points. 
The distance to the first grid point away from the wall is 10 -5 local chord lengths. The struc- 
ture of the mesh is shown in Fig. 16a. There is considerable skewness in the wing-tip region, 
as indicated in Fig. 16b. A somewhat coarser mesh has been also used with 193x49x33 
points. In the three-dimensional computations done so far, only a single time step is executed 
on the coarse meshes. Fig. 17 shows the convergence behavior of the three-dimensional 
scheme for M <*=(). 84 and Re^ll.xlO 6 . Two flow cases, namely a = 3.06° with attached flow, 
and a = 6.06° with a strong separation on the upper surface are considered. For the 
attached-flow case the influence of the fourth-difference dissipation on convergence is 
displayed in Fig. 17a. The convergence is not as good as for the two-dimensional cases. The 
largest residuals occur at the wing tip, where the grid is highly skewed. Therefore, a consid- 
erable improvement of the convergence can be expected, if the grid quality were improved in 
this region. The effect of multiple time steps on the coarse meshes for the attached-flow case 
is shown in Fig. 17b. The decrease of residuals versus multigrid cycles is not much affected 
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by multiple time stepping on the coarse meshes. Apparently, there are high frequency oscilla- 
tions on the fine mesh, which are only slowly damped due to the grid skewness, and this 
behavior is not improved using more work on the coarse meshes. The convergence of the glo- 
bal flow field, however, is improved as shown in Table 1. The improvements which are 
obtained by multiple time stepping on the coarse meshes, are more pronounced for the 
separated flow case as shown in Fig. 17c. Here, the interaction process between the shock 
and the separating boundary layers is considerably accelerated, and converged values of lift 
and drag are obtained more rapidly, see Table 1. 

Before examining details of the flow solution, it may be appropriate to look at the 
influence of the far-field distance on the lift and the drag of the wing. For this purpose, the 
distance to the far field is varied by omitting grid planes at the far field. The effect of deter- 
mining the free-stream quantities by superimposing the flow field of a horse-shoe vortex to the 
onset flow is presented in Fig. 18. For far-field distances of 7 half spans, the results with and 
without vortex differ only by 0.3%. The variation of the distance to smaller values shows, 
that the use of the vortex reduces the sensitivity of lift and drag to the distance by about 50%. 
In the results given below, the distance to the far field is kept at 7 half spans, and the horse- 
shoe vortex terms are included. 

The pressure distributions along several spanwise stations of the wing are displayed for 
the attached-flow case in Fig. 19. The results of the 289x65x49 mesh agree well with those 
from the coarser mesh and with experimental data [20]. From Table 1, it is evident, that the 
global features of the flow converge rapidly. Indeed, plots of the solution after 50 multigrid 
cycles, which are displayed in Fig. 20, show virtually no differences to the fully converged 
solution. 

Results for the separated-flow case are displayed in Fig. 21. With the B.-L. turbulence 
model, large discrepancies between computation and experiment occur. The size of the 
separated region is underpredicted, and, consequently, the shock is located too far down- 
stream. The agreement between computation and experimental data is greatly improved with 
the J.-K.-mod. model, Fig. 22. The size of the separated region and the location of the shock 
compare well with the data. The wall streamlines show a mushroom-type structure which is 
typical for wings at high angle of attack. Table 1 shows, that approximately 80 multigrid 
cycles on the fine mesh are sufficient for convergence to engineering accuracy. If the 
193x49x33 grid is considered to be sufficiently fine, 80 multigrid cycles require less than 40 
minutes CPU-time on the Cray-2 computer. 

Concluding Remarks 

A cell-vertex scheme for the Navier-Stokes equations, which is based on central difference 
approximations and Runge-Kutta time stepping, has been described. Using local time step- 
ping, implicit residual smoothing with locally varying coefficients, a multigrid method, and 
carefully controlled artificial dissipative terms, very good convergence rates are obtained for a 
wide range of two- and three-dimensional flows. Details of the acceleration techniques, which 
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are important for convergence on meshes with high-aspect-ratio cells, have been discussed. In 
general, engineering accuracy is obtained within 40 multigrid cycles on the fine mesh for 
two-dimensional flows, and within 50-80 multigrid cycles in three dimensions. The accuracy 
of the code is examined by grid refinement studies and comparison with experimental data. 
For an accurate prediction of turbulent flows with strong separations, a modified version of 
the nonequilibrium turbulence model of Johnson and King is introduced, which is well suited 
for an implementation into three-dimensional Navier-Stokes codes. It is shown that the solu- 
tions for three-dimensional flows with strong separations can be dramatically improved, when 
a nonequilibrium model of turbulence is used. 
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Table 1 : Convergence behavior for two- and three-dimensional flows. Numbers denote the 
multigrid cycles on the fine mesh required to obtain the convergence criterion. Numbers in 
parenthesis are the total computation times in seconds on Cray-2. Coefficient of the fourth 
difference dissipation K (4) = 1/32. 
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Fig. 5 . Convergence history, streamlines, and distributions of pressure and skin friction for 
laminar flow around NACA 0012 airfoil, M„.=0.5, a=0°, Re.=5000, grid 257x65. 
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Fig. 6 : Convergence behavior for transonic flow around RAE 2822 airfoil, M„=0J3, 0£=2.79° 
Re^.SxlO 6 ). 
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Fig. 7 : Influence of grid density and artificial viscosity on global force coefficients for flow 
around RAE 2822 airfoil (M„=0.73, a=2.79°, Re„=6.5xl0 6 ), computed with B.-L. model. 
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Fig. 8 : Distributions of pressure and skin friction for different grid densities (RAE 2822 air- 
foil, ^=0.73, 0=2.79°, Re„=6.5xl0 6 , K (4) =l/64), computed with B.-L. model . 
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Fig. 9 : Effect of the turbulence model on distributions of pressure and skin friction for RAE 
2822 airfoil (M„=0.73, cx=2.79°, Re oo =6.5xl0 6 , grid 481x97). 




Fig. 10 : Effect of the turbulence model on distributions of pressure and skin friction for RAE 
2822 airfoil (M m =0.75, og=2.81°, Re.=6.5xl0 6 , grid 481x97). 
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Fig. 12 : Comparison of lift coefficient versus angle of attack for the NACA 0012 airfoil, 
M.=0.7, Re^.OxlO 6 , grid 321x65. 
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Fig. 13 : Comparison of lift versus drag polars for the NACA 0012 airfoil, ^=0.7, 
Re^.OxlO 6 , grid 321x65. 
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Fig. 14 : Convergence history and distributions of pressure and skin friction for turbulent 
supersonic flow around NACA 0012 airfoil, ^=2.0, o=0°, Re^lO.OxlO 6 , grid 321x65, com- 
puted with the B.-L. model. 
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a) Influence of fourth-difference dissipation 
on convergence, single time step on coarse 
meshes, o=3.06°, B.-L. model 



b) Convergence for two time steps on 
coarse meshes, K (4) =l/32, oc^3.06°, B.-L. 
model. 



c) Convergence for nonequilibrium J.-K.-mod. turbulence model, K (4) =l/32, a^6.06°. 


Fig. 17 : Convergence behavior for transonic flow around ONERA-M6 wing M oo =0.84, 
Re^plLOxlO 6 , mesh 289x65x49 
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Fig. 18 : Lift and drag coefficients versus distance to the far-field boundary for transonic flow 
over ONERA-M6 wing M.^.84, a = 6.06°, Re^ll.OxlO 6 . 
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289x65x49 grid, 50 cycles C L = 0.2678, C D = 0.01781 
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Fig. 20 : Comparison of pressure distributions on ONERA-M6 wing, ^=0.84, a = 3.06°, 
R e oo=l l.OxlO 6 , K (4) =l/32, B.-L. model, after 50 multigrid cycles with the fully converged 
result, two time steps on the coarse meshes. 
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21 : Pressure distributions and wall streamlines for ONERA-M6 wing, M^O.84, 
6.06°, Re =11.0xl0 6 , K (4) =l/32, computed with the B.-L. model. 
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